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Abstract 

Hartree-Fock-Bogoliubov calculations of hot fission in 240 Pu have been performed with a newly- 
implemented code that uses the D1S finite-range effective interaction. The hot-scission line is 
identified in the quadrupole-octupole-moment coordinate space. Fission-fragment shapes are ex- 
tracted from the calculations. A benchmark calculation for 226 Th is obtained and compared to 
results in the literature. In addition, technical aspects of the use of HFB calculations for fission 
studies are examined in detail. In particular, the identification of scission configurations, the sen- 
sitivity of near-scission calculations to the choice of collective coordinates in the HFB iterations, 
and the formalism for the adjustment of collective-variable constraints are discussed. The power 
of the constraint-adjustment algorithm is illustrated with calculations near the critical scission 
configurations with up to seven simultaneous constraints. 
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I. INTRODUCTION 



The last three decades have seen a resurgence of interest in the microscopic description of 
nuclear fission. This renaissance in fission theory has been ushered in by progress in formal 
many-body theory and by the advent of faster and parallel computers. The microscopic 
approach can boast a well-established track record^of^accomplishment over the last three 
decades, such as t 



with temperature 
fission-isomer lifetimes 
yields 



ie prediction of fission barriers l[ 2], [], 0, (3, Q], and their evolution 
and angular momentum [2j, the prediction of fission times 



and 



the description of hot and cold fission [2[, the prediction of fission 



n 



121 ]. the description of cluster radioactivity as very asymmetric fission |l3|], and most 



recently, the calculation of fission-fragment properties (e.g., excitation energy, shape, kinetic 
energy, emitted-neutron multiplicity, angular momentum) Q, [l5]. Despite these successes 
however, the microscopic description of fission remains one of the most difficult challenges 
in nuclear physics. 

On the other hand, the promise of a microscopic theory that can reliably predict nearly 
all aspects of fission within a single, self-consistent framework is tantalizing. A fully self- 
consistent, dynamical approach to fission has been developed by the group at Bruyeres-le- 
Chatel {2, Q, Q| j and is being implemented at Livermore This approach treats both 
static and dynamic aspects of fission self-consistently and requires as its only phenomeno- 
logical input the effective interaction between the nucleons. 

A Hartree-Fock-Bogoliubov (HFB) code is the central tool for the description of the 
static aspects of fission in the microscopic method. The use of a finite-range effective inter- 



17l |. allows for the treatment of pairing within the HFB 



action, such as the D1S interaction 

formalism H in a Mly se.f— m an„er, and witaont the need for add.tiona, phe- 
nomenological parameters. The HFB calculations can be constrained by a judicious choice 
of collective variables to explore those nuclear shapes that are relevant to fission. Such 
constraints have confirmed the richness of fission phenomena, for example by revealing the 
full range of fission modes from hot (fragments formed in maximally-excited states) to cold 
(fragments formed with no excitation energy) {2]. 

In the dynamical component of the microscopic theory, a wave packet is built from 
HFB solutions constrained over all relevant nuclear shapes using the Time-Dependent 



Generator-Coordinate Method (TDGCM) 



19 



20 



21 



221 ] . In practical applications, the 
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Gaussian-Overlap Approximation (GOA) to the TDGCM can be used to produce a collec- 
tive Schrodinger equation, and therefore a collective Hamiltonian, constructed entirely from 
the single-particle degrees of freedom. The TDGCM formalism describes the nucleus in its 
lowest-energy state, as well as its collective excitations and can be extended to 



include intrinsic excitations as well 



on the way to scission. These intrinsic excitations 



are needed for a microscopic description of fission tha t g oes beyond the standard adiabatic 



approximation usually adopted in fission calculations 



26j. This comprehensive program for 



the microscopic description of induced fission has already shown the importance of dynam 



121 ] . but a great 



ical effects in the prediction of fission times [j] and fission-fragment yields 
deal of work remains to include all the relevant physics aspects in the calculation. In par- 
ticular, a detailed and quantitative understanding of scission itself remains to be developed 
even at the level of the static calculations. 

In this paper, we focus on the static aspect of the microscopic theory with three goals 
in mind: 1) to introduce the newly-developed HFB code FRANCHBRIE [la | ; which uses 
a finite-range effective interaction, 2) to examine in detail some basic technical aspects 
of fission calculations with an HFB code, and 3) to present first-time results of scission 
properties for the hot fission of 240 Pu. In section [Til we review the HFB formalism and 
discuss in detail some features of the one-center deformed harmonic-oscillator basis, formal 
and practical aspects of HFB fission calculations with multiple constraints, as well as the 
HFB convergence algorithm itself. In section IHIl we benchmark our HFB code against two- 
center calculations of scission properties for 226 Th by Dubray et al. We then apply 
the code to the identification of hot-scission configurations in 240 Pu, and the shapes of the 
nascent fragments just before scission. 

II. THEORY 

A. General HFB formalism 

For convenience, we recall the main points of the HFB formalism with a finite-range 
effective interaction in this section and refer the reader to the literature for further details 
see, e.g., 



16]. 



27 



29]). We have implemented this formalism within the code FRANCHBRIE 
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We start from the many-body Hamiltonian in second-quantized notation (see, e.g., chap- 
ter 5 in (27]), 



with the antisymmetrized two-body matrix elements 



Vmnpq = (mn \ V\ pq) - (mn \V\qp) 



and the usual anticommutation rules for particle operators 



(1) 



In thispaper, we use a finite-range effective interaction which in coordinate space takes the 



form 



28] 



V(fi,f 2 ) 
2 



J2 (Wt + B.A - H t P T - MAPt) 
»=i 

+iW L s Via x 5 (rl - r 2 ) V i2 ■ (o\ + <r 2 ) 
+t f 1 + x P CT ) 5 (n - r 2 ) f^4^ 1 + 1 ( 



2/„2 



■^)7f 



Coul 



(2) 



where V 



[2 



Vi — V 2 , V12 = Vi — V 2 , P a is the spin-exchange operator, and P T is 



the isospin-exchange operator. The Coulomb interaction Vcoui is added if both particles are 
protons, and p (f) denotes the total nuclear density. The D1S effective interaction 2, J] has 
been used for the present calculations. Given the computationally-intensive nature of the 
calculations, we have omitted contributions from the spin-orbit and Coulomb interactions to 
the pairing field. This approximation is well justified in the case of the spin-orbit interaction 
whose intensity in the singlet-even channel is very weak, but less so for the Coulomb term 
that can significantly reduce the pairing correlations for proton pairs \s0\. We note also that 
the density-dependent part of the interaction is adjusted to cancel in the singlet-even channel 
by setting Xq — 1. Consequently, only the Gaussian terms contribute to the pairing field, 
which permits the fully self-consistent application of the Bogoliubov formalism, without the 
need for arbitrary truncations of the space or the use of ad-hoc pairing forces. The Coulomb 
exchange contribution has been treated in the Slater approximation, and the two-body 
contribution to the center-of-mass correction has been included in the mean field. 
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The Bogoliubov theory jjjs| takes into account, in an approximate way, two-body corre- 
lations beyond the mean-field restriction to particle-hole excitations. The approach defines 
quasiparticle creation and destruction operators as linear combinations of the particle cre- 
ation and destruction operators, 



n 



(3) 



Assuming there exists a vacuum of the destruction operators r/^, denoted by |5), we identify 
it as the ground state of the nucleus and its energy can be written simply as a functional of 
the density matrix and the pairing tensor or, equivalently, as a functional of the generalized 
density 



R = 




R u R 12 
R 21 R 22 



(4) 



We recall that the unitarity condition of the transformation in Eq. (3) is equivalent to 



R 2 



R 



(5) 



we will therefore write the energy as 



E (p, k, A p , A n , A) 
E (p, k) - \ (5 



N, 



A„(0 



-Tr [A (R 2 - R)] 



(6) 



where E (p, k) is the expectation value of the Hamiltonian in the quasiparticle ground state, 
A p and X n are the Lagrange parameters needed to impose the appropriate average number 
of protons and neutrons, respectively, given by the matrix R. The matrix A of Lagrange 
parameters is needed to satisfy Eq. ((5|). Thus the determination of the fundamental nuclear 
state amounts to finding the generalized density matrix that minimizes Eq. j6]). Some au- 
thors recognize Eq. ((ED as the equation of a multidimensional surface, and seek its minimum 
directly using standard mathematical techniques to find the minimum of a function. Among 



these approaches, we cite the gradient method 



conjugate gradient method 



3lfl or an improved variant known as the 



32j. The number and diversity of applications using this method 



speak to its effectiveness {3, [§], Q, Q, [yj]. In our approach to the minimization of Eq. 
([6j), we start with the variational principle, 



5E (p, k, X p , A n , A) = Tr { [ft - (Ai2 + i?A - A)] <5i?} 
= 



(7) 



V5i? where 



ft 4 -' 



SE(p, k, Xp, X T 



SB 



■31 

nm 



(8) 



Taking into account Eq. 
the Bogoliubov equation 



it is possible to eliminate the constraint matrix A, leading to 



[H (R) , R] = 



(9) 



The Bogoliubov matrix Ti in Eq. ((HI) is constructed with the help of the block matrices 
defined by Eq. (jSJ). The explicit form of these matrix elements for the D1S effective inter- 



action is given by references 



28 



The solution of Eq. (J9J) is then found by successive 
diagonalizations of the Bogoliubov Hamiltonian. This iterative solution method is described 
in greater detail in section HTDl and appendix lAl 



B. Basis truncation and aspects of one-center basis calculations 



In practical applications, the formalism of section III Al must be expressed in some basis. 
Typically, the deformed Harmonic-oscillator (HO) basis (see, e.g., chapter 2 in [27|) has been 
used in many HFB calculations, including those dealing with fission |4|, |36|. The basis states 
in cylindrical coordinates (p, z, <p) are 



(f\n r ,A,n Z} a) = $„ r ,|A| (p; b±) - 
x$„ z (z; b z )xc 



/27T 



(10) 



where the explicit forms used in this work for the radial (<Ey.,|A|) and Cartesian (<& n J com- 



ponents and their relevant properties can be found, e.g., in [35|, and x<? is a spinor function 
for a = ±1/2. These basis states assume axial symmetry of the nucleus explicitly. Other 
symmetries can also be imposed on the HFB calculation to reduce the overall size of the 
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problem. Two symmetries in particular are relevant to the fission calculations in this paper: 
the symmetry with respect to the parity operator II 

TL\n r ,A,n z ,a) = (— l)' A ' +nz \n r , A,n Zl a) 

and the symmetry with respect to the z-signature operator S z = iR z (tt), where R z (ir) effects 
a rotation by it in both spatial and spin space, 



S z \n r , A,n z ,a) = cr(-l) |A| \n r ,A,n z ,a) 

Throughout this work, only the z-signature symmetry has been imposed, leaving the fis- 
sioning nucleus free to violate the symmetry with respect to parity and assume asymmetric 
shapes. These symmetries are taken into account explicitly by rewriting the general Bogoli- 
ubov transformation of Eq. (J3]) in terms of the relevant quantum numbers as 

vl( q ,s z ,Q) = J2[u^l(q,s z ,Q) 



+V q > Sz 

' nfx 



(q, s z , tt) 



E 

n 



u q 



a n (q,s z ,Q)] 
" <A a n (q,s z ,n) 



+ ( Vtf" 



4 (?,«*, fi) 



where q distinguishes protons and neutrons, s z = ±1 is the z-signature quantum number, 
and fl is the total angular-momentum projection for the time-reversed state. 

Even with the z-signature symmetry imposed, the treatment of fission can require large 
basis sizes and the calculation of a large number of two-body matrix elements. In order 
to further limit the size of problem, various basis truncation schemes have been devised. 



Some [37|| keep only those basis states with corresponding HO energies below a given cutoff, 
while other schemes 4, 3^| directly allow for more quanta along the z direction-the direction 
of elongation of the fissioning nucleus-compared to the radial direction. In the truncation 



scheme of 



371 ], the HO quantum numbers must satisfy 

hw L {n± + 1) + hio z (n z + - ] < hu (N + 2) 



(11) 



with n± = 2n r + \A\ and for a given maximum shell number N, where the oscillator frequen- 
cies are related to the length parameters b± and b z in Eq. (TTOl) by 



h 

mb 



2 ■ 



h 

mbl 



WiWz 



(12) 



and m is the nucleoli mass. With increasing axial elongation and for fixed N, Eq. (PTTT) adds 
more shells in the z direction while simultaneously decreasing the number of shells in the 
radial direction, thus keeping the basis size from growing too quickly with deformation. In 



the truncation scheme of 



4 



381 ] . the condition 



n 



— + 2n r + |A| < N (13) 

q 

is imposed for a given maximum shell number N and parameter q. In this work we have 
used both truncation schemes. The truncation given by Eq. ([Til has been used for most 
calculations in this paper, while the truncation of Eq. ( fl~3l ) has been used mainly in section 

una 

The oscillator lengths b± and b z in Eq. (fTOj) . or equivalently the frequencies u± and u z , 
are variational parameters in the HFB calculation that must be chosen to minimize the HFB 
energy. Through a series of calculations in 240 Pu using the truncation scheme of Eq. ffl~3l) 
with N = 13 and q = 1.5, and exploring a wide range of values of the constraints on the 
quadrupole (Q20) and octupole (Q30) moments, an approximate dependence was obtained 
for the frequencies that minimize the HFB energy, given by 

fko = 8.4345 - 0.0021668 Q 20 (14) 

— = 1.7041 + 0.0028743 Q20 (15) 

u g 

with Q20 in barns and ftw in MeV. No significant dependence on Q 30 was observed in the 
range of interest. 

Perhaps the most important aspect of the basis states in Eq. ( fTOl) is that they are centered 
about the origin by construction. In particular, the Gaussian factor in Eq. (fTOl ) ensures that 
the nuclear wave function falls off rapidly with increasing z. Despite this feature of the basis 
states, we will show that it is still possible to describe the exotic shapes occurring in fission. 
In order to describe both the neck (near z — 0) and nascent fragments (typically 5-10 fm 
from the origin) with the basis states of Eq. (fTOl) . we are forced to include many quanta in 
the z direction, and to use relatively large values of b z . 

To justify the use of the one-center basis for the range of fissioning configurations and 
quantities examined in this paper, we have performed separate HFB calculations for 134 Te 
and 106 Mo centered at the origin, and translated the resulting wave functions to the typical 
positions these nuclei occupy as Pu nascent fission fragments. The formalism required for 
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Figure 1: (Color online) Plots of the nuclear densities for fragments of 134 Te and 106 Mo along the axis 
of elongation of the nucleus, calculated in the one-center basis and plotted (as solid black lines) cen- 
tered at z = -7.63 and 9.65 fm, respectively. The dashed red lines represent the same densities, but 
translated from the origin to their respective centroid positions within a finite harmonic-oscillator- 
basis truncated according to Eq. (fTTj) and with iV=13 shells, using the formalism in appendix 

m 

translating a wave function expressed within a finite HO basis is given in appendix [Bl The 
basis was truncated according to Eq. ffTTjl with N = 13, resulting in a maximum number 
n z = 26 along the z axis. The result is shown in Fig. [Q, and compared to a translation in an 
infinite-sized basis (obtained in practice by redrawing the curves at the displaced centroid 
positions while preserving their shape). The comparison clearly shows the appearance of 
spurious tails for each fragment translated within a finite-size basis. If the fragments are 
separated further, e.g. by an additional 2.5 fm for each fragment in Fig. O the tails grow 
larger. However, the tails caused by the translation in a finite basis remain relatively small 
(~ 10~ 4 fm -1 in Fig. [Q, and ~ 5 x 1CT 4 fm -1 in Fig. [2]), and the separations between the 
fragments in both figures are larger than those encountered in the remainder of this work. 
In section [HI Al we will show that these tails do not significantly affect the nuclear properties 
calculated in this paper. In a forthcoming publication jsjj] we will explore a more microscopic 
definition of scission and of the fission fragments, and we will calculate quantities such as 
the interaction energy between the fragments that may be more sensitively affected by the 



presence of these tails 



43|]. 
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Figure 2: (Color online) Same as Fig. Q], but for the 134 Te and 106 Mo fragments translated an 
additional 2.5 fm each, to centroids at z = -10.13 and 12.15 fm, respectively. 

C. Multiple constraints in HFB calculations 

In this section, we focus on formal and practical considerations in the choice and control of 
multiple constraints in HFB calculations. We will describe a mechanism for the adjustment 



of the constraints which generalizes the discussion in [28fl. The formalism described here and 
used in our calculations is that of variation with linear constraints. Other approaches for 
the adjustment of constraints, such as the quadratic-constraint method can also be found 
in the literature {40}. We have adopted the linear- variation approach in our work because 
we have found it to be stable and robust, and these are important qualities needed to map 
out the scission configurations, which requires precise control of the nuclear shape. For a 
process like fission, these constraints are central not only to being able to drive the nucleus 
to scission, but also to uncover the full richness of the microscopic method in its ability to 
describe the complexities of fission. In section III Al we already discussed the introduction 
of constraints on the average number of neutrons and protons for the HFB Hamiltonian. 
Further constraints can be introduced through the external-field one-body operators AjFj, 

H-J2*A (16) 



10 



where the parameters Aj are used to adjust the field intensities. Based on Eq. j9]), the 
Bogoliubov equation associated with Eq. (flB)) can now be written 



H(R)-J2 A ^ R 



where 




(17) 



in the particle-hole representation, and TC (R) is given by Eq. (JSj) . In what follows, we will 
use the notation 

H(R,{\ t }) = H(R)-Y,^i 

i 

where {Aj} represents the set of Lagrange multipliers other than those associated with the 
proton and neutron numbers. The A« Lagrange multipliers can be adjusted to yield an HFB 
solution with desired expectation values fi of the fields 

(fy = ^TrFj + ^TrWiR 

= fi 

The formalism used to find the appropriate A, parameters is derived in appendix [A] In 
describing fission within the microscopic approach, we are free to impose any number of 
constraints, each defined by a corresponding external-field operator. We are limited in this 
task by the computational requirements, which grow quickly with the number of constraints, 
and by their relevance to the fission process. 

In the simplest physical picture of fission, we expect that the nucleus will stretch along its 
symmetry axis until scission, and therefore introduce the mass quadrupole operator Q20 as 
a constraint. Next, the octupole operator Q30 is introduced to account for the range of mass 
divisions observed in fragments, from symmetric to asymmetric. With the introduction of 
the octupole constraint, we are forced to impose a constraint on the dipole moment, Qio, as 
well in order to maintain the center of mass of the nucleus fixed. The hexadecapole operator 
Q40 controls the formation of the neck between nascent fragments, and accounts for the 



range of fission modes from cold to hot {2]. In addition, we recall that the HFB procedure 
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Figure 3: Calculated HFB energy for 240 Pu as a function of hexadecapole moment, and for 
quadrupole moments of 300 b (cold fission) and 370 b (hot fission). For the (Q20) = 300b case, the 
fission valley is seen near (Q4o\ = 130b 2 , and the fusion valley is near ^Q 40 ) = 90 b 2 . For the 
(Q20) = 370b case, only the fusion valley is observed, near (Q40) = 140b 2 . 



requires constraints on the expected values of the proton-number (N p ) and neutron-number 
(N n ) operators. 

In Fig. [3l we show a calculation of the HFB energy for 240 Pu as a function of Q40 
(Q40 = (q^) at two quadrupole deformations, 300 b and 370 b, which correspond to the 
so-called cold and hot fission limits, respectively [2]]. These calculations were performed with 



5 constraints (for the values of (Nj 



94, (N, 



146, (Q 



'10 



0, (Q 



{20 



300 b or 



370b , and 80 b 2 < {Q m J < 200 b 2 ). In the cold-fission case, a barrier of height ~ 4.0 MeV 
relative to the fission-valley minimum separates the two valleys. Near the hot-fission limit, 
the fission valley has disappeared and the nucleus spontaneously falls into the fusion valley 
near (^Q^ = 140 b 2 . Between the hot and cold extremes, the nucleus can undergo fission 
through a range of intermediate modes. 

The energy curves plotted in Fig. [3] effectively represent slices at fixed values of (^Q^ 
in Fig. 3 of [2j. The most striking feature in Fig. [3] is the sudden variation in energy over a 
very small step size in (Qio^ of 1 b 2 . In the cold-fission case, a drop of 2.7 MeV is observed 
in going from (^Q^ — HOb 2 to 109 b 2 , and in the hot-fission case a more pronounced 
drop of 7.6 MeV occurs in going from = 190 b 2 to 189 b 2 . These abrupt changes in 

energy, which are in contrast to the smooth behavior displayed in j^, correspond to a sudden 
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Figure 4: (Color online) Calculated nuclear densities in steps of A (^Q^oJ = 1 b 2 around the scission 
configuration for cold (top panel) and hot (bottom panel) fission. The legends give the values of 
Q40) for the different curves. 



reduction in the neck size (Fig. 0j), which we take as an indicator of a transitional phase 
where the nucleus is undergoing scission. Note that the identification of such transitional 
phases requires extremely small variations of the constraints, which could explain why they 
were not seen in |2|. The precise control of the constraints needed to study the region around 
scission is one of the important points that emerges from the work presented in this paper, 
and the motivation for going into some detail in the description of the constraint-adjustment 
algorithm in the next section and in appendix lAl 

The rapid change of the neck size mentioned above suggests the introduction of a con- 
straint proportional to the average number of particles (QnJ in the neck separating the 
nascent fragments, where jj| 

\2" 



Q N = exp 



(z - z N y 



(18) 



with a^v = 1 fm, and z^ is the position of the neck (defined as the point between the 
fragments where the matter density is lowest). As shown in Fig. the energy calculated as 
a function of \ QnJ becomes smoother and continuous. A more detailed discussion of this 
result is given in the latter part of section III Dl 
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Figure 5: Variation of the HFB energy as a function of the number of particles in the neck, 
defined by Eq. (fl~8l) . at the scission configuration ((^Q^o^ = 189b 2 ) for the hot-fission calculation 
.Q20) = 370b) in Fig. 1 



D. The HFB convergence algorithm 

The control of HFB calculations with multiple constraints is a delicate procedure, made 
difficult by the number of constraints and their inherent correlations. Because the topic 
continues to be of current interest in problems that rely on constrained-HFB methods even 



beyond fission 



32 



411 ] . the convergence algorithm used in the present HFB calculations is 



discussed in detail here. The algorithm must balance, at each iteration, the diagonalization of 
the HFB Hamiltonian to ensure self-consistency, and adjustment of the Lagrange multipliers 
in Eq. (fi~6l) . The main steps of the algorithm are as follows 



1. Read initial generalized density R and Lagrange multipliers Aj 

2. Construct constrained HFB Hamiltonian H (R, 

3. Diagonalize Ti (R, {Xi}) 



4. Construct new R 



5. Mix R between consecutive iterations using a mixing parameter a (see Eq. (1201) ) 



6. Adjust value of a based on convergence criterion 
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7. 



8. 



9. 



Calculate SXi needed to yield desired constraint values, adjust Aj 
Calculate 5R corresponding to the SXi, adjust R 
If HFB solution is not converged, return to step [2] 



The first 4 steps in this algorithm are fairly self-explanatory and make use of the formalism 
derived in section III AL We will examine the remaining steps in greater detail since they are 
not typically discussed in depth in the literature. 

At the end of each iteration i, the convergence of the HFB solution is assessed by cal- 
culating the largest variation from the previous iteration in the elements of the generalized 
density matrix, 

e, = sup\R% n (i)-R% n (i-l)\ (19) 

The quantity Si is also used to determine the coefficient a in step [5] which mixes the gener- 
alized densities between successive iterations using an adjustable coefficient a, 

iC n (0 - (1 - a) R P Z (<) + aiC (< - 1) (20) 

with < a < 1. This mixing is essential to slow down the convergence algorithm which 
would otherwise often behave erratically in the first few iterations and could fail to converge 
at all. The mixing coefficient a is adjusted in step [6] in such a way that it tends to zero as £j 
decreases. In practice, two thresholds are supplied, e m i n and e max , along with a maximum 
value a max for the mixing coefficient such that 

Oi = < ry £ i- £ min - . p. ^ p 

i ^max _ _ . °min ^ °i ^ °max 

tmax tmin 

Furthermore, if the HFB solution diverges from one iteration to the next (i.e., if > £i_i) 
then a is set to a max and remains at that value until the HFB solution converges again. For 
the work in this paper we have used e m i n = 1CT 3 or 10~ 4 , e max = 10 _1 , and a max = 0.5 
(or in a few cases 0.8 for a slower initial convergence). We note in passing that the mixing 
of generalized density matrices is a global operation, i.e. the same coefficient a is used 
for all the matrix elements. The Broyden method, or its more elaborate modified version 
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Figure 6: Plot of the convergence metric, given by Eq. (TT9|> . as a function of HFB iteration number 
for the ((340/ = 110 b 2 cold-fission point in Fig. [3l 



411 ]. could provide a better alternative for optimizing the choice of the mixing coefficient by 
associating an independent value of a to each matrix element. 

The formalism needed to adjust the Lagrange parameters in step [3, and the generalized 
density in step [8] is presented in appendix El and we stress the importance of adjusting both 
for a stable convergence of the HFB method. The algorithm is considered to have converged 
in step M if Si < £ m in for several iterations (typically 2 in the present work). 

In order to illustrate various aspects of the convergence algorithm, we have examined the 
cold-fission point at (Qao^ = 110 b 2 in Fig. [3] in detail. Because this point corresponds 
to a local maximum in the HFB energy, its calculation is particularly demanding on the 
convergence algorithm. In Fig. [6] we show the convergence criterion, e, calculated using 
Eq. (fT9l) at each iteration. The HFB solution is found to better than e < 10~ 4 after 156 
iterations in this case. We note a region in Fig. [6] roughly between iterations 10 and 40, 
where e appears to be relatively constant and the convergence is correspondingly slow. In 
this region, all the constraints appear to be close to their desired values, except for the dipole 
moment. The (QiqJ value is still relatively large (~ 0.06 — 0.2 fm) and may be responsible 
for the stagnant convergence. 

In Fig. [7] we examine the adjustment of the five constraints at each iteration. The 
figure shows the relative deviation of each constraint from the desired value. For all but the 
dipole-moment constraint, this relative deviation of the calculated average value (Q) of the 
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Figure 7: (Color online) Relative deviations of the calculated constraint values from their desired 
values as a function of HFB iteration number for the calculation with (^Q^ = 110b 2 . The relative 
deviation for the dipole moment is given by Eq. (|22l) . and by Eq. (|21~j) for all other constraints. 
The constraints shown are: (Qio^ (black solid line), (Q20^ (red dotted line), (^Qio^ (green dashed 
line) , ( N n \ (blue dot-dashed line) , and ( N p \ (turquoise dot-dot-dashed line) . 



constraint from its desired value q is given by 

'Q)-q 



q 



(21) 



In the case of the dipole moment, the desired value is gio = and Eq. ( 1211 1 cannot be 
used. Instead, we obtain from (QioJ the position of the centroid of the nucleus, given by 
(Qio) I A- where A = 240 is the total number of nucleons, and compare it to the calculated 
root-mean-squared radius of the nucleus, R rms , using the ratio 



Q 



10 



AR r 



(22) 



The calculation is started from an HFB solution that differs only in the value of the 
hexadecapole constraint, (Qio^ = 115 b 2 , with all other constraints the same. Hence we see 
in Fig. [7]that at the first iteration, all relative deviations except the one for the hexadecapole- 
moment constraint are small. The calculation converges to the desired level of accuracy after 
156 iterations. 
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Figure 8: (Color online) Same as Fig. [7j but for the calculation with (Q40) = 130 b 



This difficult convergence should be contrasted with the calculation of the cold-fission 
point at = 130 b 2 , near the bottom of the fission valley in Fig. [3j The relative 

deviations of the constraints for this more stable calculation are shown in Fig. El After 
the tenth iteration, all constraints tend to the desired value rapidly and smoothly. This 
calculation is converged to the same level of accuracy as the one at (Qm^ = 110 b 2 after 
only 33 iterations. 

Finally, we discuss in greater detail the discontinuities observed in Fig. [3l Such discon- 



tinuities have been alluded to in the literature |42j as a potential difficulty for microscopic 
calculations. In this section, we show how these discontinuities are an indicator of a change 
in the meaning of certain collective coordinates near the critical scission configurations. We 
also show how these discontinuities can be eliminated through the choice of a more appro- 
priate collective coordinate. 

The impact of these discontinuities can be felt even before the scission configuration is 
reached. We illustrate this point by showing the results of HFB calculations, performed 
with identical multipole constraints up to the hexadecapole moment (i.e., with the same 
(Qio^i values), but different initial densities. We will approach the 

cold-fission scission configuration near (Q40) = HOb 2 in Fig. [3] with an initial density 



corresponding to either a scissioned or non-scissioned nucleus. The first calculation, shown 
in Fig. [H was performed at (Q40 / = 130 b 2 , near the bottom of the fission valley. Two curves 
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are shown, corresponding to a initial choice of the generalized density calculated at = 
135b 2 (whole nucleus), and = 90b 2 (broken/scissioned nucleus). As expected, both 

choices of starting point lead to exactly the same HFB solution, as is evidenced by the 
overlapping density curves in Fig. [9l By contrast, Fig. [TO] compares calculations at (Q40 

,2 /• „ „„ • •„ \ „,„,• „ r i , / A \ , ni2 



115 b (i.e., near scission), starting from solutions at {^QioJ = 120 b (whole) and 
90b 2 (broken). Both solutions have the same values of the first four moments, yet the 
calculation started from a whole solution leads to a whole result, while the broken starting 
configuration leads to a broken-nucleus solution. A similar effect is observed in Fig. ITTl 
corresponding to a calculation very close to scission at = 110 b 2 with starting densities 

from — H5b 2 (whole) and (Qao^ = 90 b 2 (broken) solution. Note that these HFB 

calculations are performed with an unprecedented 7 simultaneous constraints. 

The densities plotted in Figs. IMTTl reveal a complex relationship between the hexade- 
capole and Qn degrees of freedom. These two coordinates are not related by a one-to-one 
mapping and cannot be used interchangeably to drive the system to scission. In Fig. IT2l we 
show the HFB energy surface as a function of Q 40 and Qn for the calculation with all mo- 
ments up to hexadecapole fixed. In particular, ((^20^ — 300 b, and (Qso^ — 34.951 b 3//2 -the 
value of the octupole moment for the two calculations in Fig. [TTj The shape of the en- 
ergy surface suggests that energy-minimizing HFB solutions can exist which have the same 
value of (Qao^, but distinct values of (^Qn^- For most-but not all-values of a small 

barrier in the surface (marked by a solid line along the surface in the figure) separates the 
minima with differing values of (Qn^- This barrier is at best a few hundred keV's high 
and decreases rapidly with decreasing as we approach the scission configuration. At 

= 110 b 2 the barrier has dropped to only 1.8 keV and vanishes completely between 
= 104 b 2 and 110 b 2 . This break in the barrier causes the discontinuity in Fig. [3l 
where the calculations are performed without a constraint on (Qn^ to prevent the HFB 
calculation from falling into the scissioned configuration. 

Near scission, the total multipole moments of the nucleus are determined by the in- 
trinsic and relative moments of the fragments, and rearrangements between these terms can 
produce different matter distributions with the same overall moments, at least up to the hex- 
adecapole. Thus imposing a constraint on ^40^) will not necessarily result in a constraint 
on the neck size near scission. The (Qn) constraint on the other hand was already shown 
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Figure 9: (Color online) Comparison of nuclear densities for the = 130b 2 cold-fission point 

in Fig. [3j starting either from a whole (solid black line) or scissioned/broken (dashed red line) initial 
configuration of the nuclear density in the HFB iterations. All moments up to the hexadecapole 
have been constrained to the same values for the two calculations. 
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Figure 10: (Color online) Same as Fig. El but for a calculation at /<54o^ = 115 b 2 . 

to produce a smooth energy dependence in Fig. [5] and is therefore the suitable coordinate 
in the study of fission for configurations near and beyond scission. 
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Figure 12: Energy surface calculated with constraints on y^ny = 146, \Np} = 94, yQioy = 0, 
^Q 20 ^> = 300b, (q 30 \ = 34.951b 3/2 , 90b 3/2 < ^Q 4 o) < 130b 3/2 , and 0.05 < ^Q N ^ < 3.05. The 
dark lines along the surface mark the position of a small local barrier on the surface. 

E. Scission in the constrained-HFB approach 

In this section, we briefly discuss various signatures of scission. Some of the characteristics 
of scission have already been mentioned in sections Til CI and Til Dl The standard indicators of 
scission are sudden changes in either energy (interaction energy between fragments or total 
HFB energy) or shape (neck size or hexadecapole moment) for the nucleus fl^ |. For the 
work in this paper, we use the same semiclassical definition of the nascent fission fragments 
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as in [14]], where a position along the symmetry axis of the nucleus is identified as a divider 
between left and right fragments, and the fragment properties are obtained as integrals over 
;he density with this cut as an endpoint for the integrals. In a forthcoming publication 
39], we will adopt a more microscopic criterion to identify the fragment {43]], based on the 



individual single-particle wave functions, and using the changes in the interaction energy 
between fragments as an indicator of scission. In this paper we will focus instead on the 
HFB energy and the number of particles in the neck before and after scission. 

Consider, for example, the cold-fission calculation in Fig. [3] At — HOb 2 there 

is still a significant amount of matter in the neck connecting the nascent fragment with 
(Qn) = 2.41. At (Qao^ = 109b 2 however, the neck breaks and (Qn/ drops to 0.50 
particles. This sudden variation in shape over a small increment in hexadecapole moment is 
shown in the top panel of Fig. [H At = 90 b 2 , the bottom of the fusion valley, \ Qn\ 

has been reduced to 0.09 particles. From — HOb 2 to 109b 2 , the total HFB energy 

drops by 2.7 MeV, and the difference in energy between ^Q 40 ) = 110 b 2 and 90b 2 is 10.2 
MeV. 

A similar analysis can be performed for the hot-fission calculation in Fig. [3] In this 
case, the last point where a sizable neck still exists between the nascent fragment is at 
(Q±o^ = 190b 2 , with (Qn) = 2.92 particles. By = 189b 2 the neck has essentially 

disappeared, and (Qn^ has dropped to 0.23 particles. At the bottom of the fusion valley, 
where (^Q^ = 140 b 2 , there are only (^Qn^ = 0.02 particles in the neck. The change in 
shape is plotted in the bottom panel of Fig. [H The drops in energy are more significant 
than in the cold-fission case. From ^Q 4 o) = 190 b 2 to 189 b 2 , the total HFB energy drops 
by 7.6 MeV, and from (q 40 \ = 190 b 2 to 140 b 2 , it drops by 20.1 MeV. 



III. RESULTS 



A. Benchmark: 226 Th scission 

We have performed HFB calculations of hot-fission properties for 226 Th, in order to 
compare with the results in jl^| that were obtained with two-center HFB calculations. We 
have used both the basis truncation of Eq. (fT3l with N = 13 and q = 1.5, and the one 
given by Eq. (fTTj) with N = 13. The oscillator-frequency parametrization of Eqs. (fT4l) and 
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(fT5l) was used, even though it was obtained for calculations in 240 Pu. We will show that our 
results are in good agreement with those of Dubray et al. for 226 Th with either basis 
truncation scheme. 

In Fig. [CH we plot the hot-scission line for 226 Th, and compare it to the one obtained in 
\\\ . The scission line was determined by performing series of calculations at fixed 
and increasing values of (Q2o^ by 5 b, each calculation using the previous one as a starting 
point, until an HFB solution was found where the neck size decreased drastically. Lines 
separated by A (Q2o^ = 5 b connecting the HFB solutions just before and just after the 
breaking of the neck are displayed in Fig. |T3l bracketing the actual scission line. These lines 
are in good agreement with the 226 Th scission line in [141 ] . In Fig. [14] we examine the region 
with (Q-so^ = 25 — 35 b 3 ^ 2 in greater detail. A series of HFB calculations were performed at 
constant (^20^ values of 280, 310, 360, and 400 b starting from (^30^) = 25 b 3//2 in each case 
and proceeding in steps of A ((^30^ = 1 b 3//2 . For these calculations, the basis truncation 
of Eq. (fTTI) was used with iV = 13 in order to provide a larger number of oscillator shells 
(up to 26 in practice) in the z direction, while keeping the overall number of basis states 
relatively low. With these large-basis calculations, we find that the results of Dubray et al. 
[14I ] are very well reproduced. 

In Fig. [T5J, we compare the mass quadrupole moment calculated for the fragments for 
the HFB solutions just before scission (solid disks connected by solid lines in Fig. [13]) to the 



corresponding result in 



141 ] . As in [14j, the Q20 values were calculated by integration over 



the left- and right-fragment densities, truncated at the neck position. The results of [14| are 
well reproduced by our calculations. Similarly, in Fig. [16], we show the octupole moment 
of the fragments compared to the Dubray et al. results. In this case as well, the agreement 
between the two sets of calculations is good. 

The agreement between one-center and two-center calculations in Figs. (fT3 l - ([T6l is reas- 
suring, both as a benchmark for the HFB code used in this work, and as an assessment of 
the applicability of the one-center basis near scission. With these results in mind, we turn 
next to the fission properties of 240 Pu. 
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Figure 13: (Color online) Scission line for 226 Th obtained in this work, and compared to the result 
of Dubray et al. |l4j|. The solid disks connected by a solid green line represent HFB solutions just 
before scission in this work, and the solid disks connected by a dashed red line represent solutions 




after scission in this work. The thick solid black curve is the scission line taken from 
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Figure 14: (Color online) Large-basis HFB calculations in 226 Th along lines with fixed ( Q 



!20 



performed to reproduce the details of the scission line found in Dubray et al. HJ|. A dashed line 
connects the last point before scission, and should be compared to the Dubray et al. result (solid 
line). 
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Figure 15: (Color online) Comparison of fission-fragment quadrupole moments as a function of 
fragment mass number between this work (solid black disks) and the results in (solid red 
triangles). 
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Figure 16: (Color online) Same as Fig. HH but for the fission-fragment octupole moments. 



B. 240 Pu scission 



For the 240 Pu calculations, we have used the truncation scheme of Eq. (TTTT) with N = 13. 
The parameterization in Eqs. (TT4l) and (TT5l) was adopted for the HO frequencies. 

Fig. [T7| illustrates the search for the hot-scission line in 240 Pu. Points along lines with 
fixed (Q-so^ or (Q2o^ increasing in steps of 1 b 3//2 and 5 b near the scission line, respectively, 
denote individual HFB calculations, each using the previous one as a starting point. As in 
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Figure 17: (Color online) Scission line for 240 Pu obtained in this work. All calculations were done 
using the basis truncation of Eq. (fTT) . The solid green disks represent HFB calculations producing 
a whole (non-scissioned) nuclear density. The empty red circles connected by a solid line represent 
scissioned configurations. 

the case of 226 Th in Fig. \T3^ the nucleus tends to stretch to much larger deformations in the 
symmetric limit. This leads to fragments that are formed much further apart in symmetric 
fission, and a corresponding drop in their mutual 



Coulomb repulsion-and therefore their 
total kinetic energy-as observed experimentally [4j|. As in the case of 226 Th, we also observe 
regions around Q 20 = 550 b/Q 30 = 35 b 3/2 and Q 20 = 400 b/Q 30 = 38 b 3/2 where the scission 
line "bulges out". In these regions, for a given Q 3Q value, the nucleus may scission at more 
than one value of Q 2 o. 

Fig. [18] compares the total HFB energy of the fissioning nucleus just before and just after 
scission. In general, scission is accompanied by a marked drop in HFB energy. That drop, 
however, is much more pronounced for fission near the symmetric limit, where it can be as 
large as ~ 50 MeV over the A (q 2 o^ = 5b change in quadrupole moment. Note that the 
fragment masses in Fig. [18] are not the same before and after scission. This difference is an 
indication of the drastic variations in the nuclear density, and the redistribution of particles 
in the neck between the two fragments at scission. 

The number of particles in the neck just before and after scission is shown in Fig. [19] as 
a function of the heavy-fragment mass. The variation in (QnJ is quite large (typically by 
an order of magnitude, but near the symmetric limit, by more than a factor of 1000). 



As in 



14 ]. we extract the fragment properties for each mass division from the HFB calcu- 
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Figure 18: (Color online) HFB energy of the fissioning nucleus, plotted as a function of the heavy- 
fragment mass number, obtained from the HFB calculations just before (solid green disks) and just 
after (empty red circles) scission in Fig. [TT1 
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Figure 19: (Color online) Number of particles in the neck of the fissioning nucleus, plotted as a 
function of the heavy-fragment mass number, obtained from the HFB calculations just before (solid 
green disks) and just after (empty red circles) scission in Fig. [TTl 



lation just before scission. However, we go further than the calculation in [IJ] by attempting 
to approach the scission configuration even more closely. We introduce an additional con- 
straint on Qn to each point in the Q20 — Q30 map of Fig. [17] just before the scission line, 
and search for the Qn value marking a point just before a drop in Ehfb occurs. Fig. EO 
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Figure 20: (Color online) Identification of the last configuration before scission for HFB calculations 
at fixed Q 30 = 10b 3/2 and 55b 3/2 , as a function of the Qn constraint. The circled points on each 
curve were chosen as the last pre-scission configuration, before the drop in HFB energy as a function 
of decreasing Qn. 



shows some typical choices for this point. In Fig. EH the charge and mass of each fragment 
is plotted, covering a range from A = 93 to 147. We note that there is a nearly linear 
relationship between the mass and charge of the fragments, which can be fitted as 

Z = 3.5349 + 0.36221 A 

This result is consistent with the prediction of the Unchanged-Charge Division (UCD) model 



45|, also shown in Fig. [21] for comparison, which for 240 Pu yields 

94 



Z 



240' 



-A « 0.3917A 



The moments of the fragments are shown in Figs. [2211241 The overall shape of the 
quadrupole moment in Fig. [22] is similar to the one shown for 226 Th in Fig. HH with a 
maximum at the symmetric limit, and a drop-off on either side. There is also a significant 
dip in the ((^20^ value near the nearly-spherical 134 Te fragment. The fragment octupole 
moment, plotted in Fig. [23], also shows similarities in shape as well as magnitude to the 
226 Th case in Fig. [16] 46j. Finally, we also show the hexadecapole moment of the fragments 
in Fig. [24] There as well, the value of reaches a maximum near the symmetric limit, 

and drops off on either side. In all cases, a line has been drawn to guide the eye using a 
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Figure 21: (Color online) Fission-fragment charge number plotted as a function of mass number, 
obtained from the HFB calculations immediately prior to scission in Fig. \T7\ The UCD prediction 
(solid red line) is plotted for comparison. 



polynomial fit to the points. The HFB calculations in Figs. 12211241 exhibit a great deal of 
fluctuation about the smooth polynomial fit. These fluctuations are due for the most part to 
the difficulty in identifying a scission configuration based on the criterion of sudden changes 



in global nuclear properties, such as the total energy. In a forthcoming paper 39], we will 
embark on a more detailed study of the scission configurations at the microscopic level, and 
extract the excitation, kinetic, and interaction energies of the fragments. The merits and 
difficulties of a scission criterion based on the interaction energy between the fragments will 
be discussed in detail. 



IV. CONCLUSION 

We have developed the HFB code FRANCHBRIE for microscopic fission studies using the 
finite-range D1S effective interaction. The code allows for the multiple constraints needed to 
explore the nuclear densities relevant to fission, and is based on matrix elements calculated 
in a one-center deformed harmonic-oscillator basis. We have provided a detailed derivation 
of the formalism required for the adjustment of those multiple constraints. 

We have applied the code to the calculation of scission configurations in the hot fission 
of 240 Pu. These calculations are relevant to studies of thermal neutron-induced fission on a 
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Figure 22: (Color online) Fission-fragment quadrupole moments, plotted as a function of fragment 
mass number, obtained from the HFB calculations immediately prior to scission in Fig. [TTl A line 
has been drawn through the HFB results to guide the eye. 
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Figure 23: (Color online) Same as Fig. [22j but for the fission-fragment octupole moments. 

target of 239 Pu. We have focused on the technical aspects of using the HFB formalism for 
fission studies. In particular, we have discussed some aspects of fission calculations within 
a one-center basis, and the importance the choice of collective coordinates in the HFB 
iterations for nearly-scissioned configurations. A scission line in the quadrupole-octupole 
plane was obtained and shows a tendency for the nucleus to reach much larger elongations 
in the symmetric limit before scission occurs. A similar feature was observed in the scission 
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Figure 24: (Color online) Same as Fig. [22], but for the fission-fragment hexadecapole moments. 

line of 226 Th by Dubray et al. [lj| using two-center HFB calculations, reproduced in this 
work with a one-center calculation. The increased "malleability" of the nucleus near the 
symmetric limit is reflected in the various moments (quadrupole, octupole, hexadecapole) 
calculated for the fission fragments and presented here. 

In a forthcoming publication, we will extract the excitation and kinetic energies of the 
fission fragments. We will introduce a microscopic criterion for the identification of fission 
fragments, and calculate their interaction energies, with special attention to the density tails 
discussed in this paper. Finally, the static calculations of hot fission presented here are the 
first step in a fully dynamical calculation of 240 Pu fission. Further developments are planned 
to explore all fission modes, from hot to cold, and to include the dynamical aspects of the 
theory in the calculations. 
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Appendix A: MULTIPLE CONSTRAINT FORMALISM 



1. Effect of the variation of a single Lagrange multiplier on the generalized density 



In this appendix, we derive the formalism for solving the HFB equation with multiple 



281 ] to the case of multiple con- 



constraints. The derivation generalizes the discussion in 
straints. 

In the first section, we give the essential formulas used in the adjustment of constraints. 
A second section illustrates the formalism with the special case of a single constraint, and 
the last section presents the general case of multiple constraints. Starting from the HFB 
equation, Eq. Q, we write for a Hamiltonian with a single constraint XF introduced as is 
Eq. (US), 



[H(R(\),\),R(\)] = 



where 



A 



A 



/(A) 

-TrF + -TrFi? (A) 

2 2 v ' 



(Al) 



is the expectation value of F in the corresponding HFB solution |A), with F given by Eq. 
( 1T71) . Consider a small variation 5X of the Lagrange multiplier, leading to a new HFB solution 
with 



[H (R (A + 5X) , A + SX) ,R(X + SX)} = 



(A2) 



where 



(X + 5X F X + 5X^ = f{X + 5X) 



= -TrF 
2 

+^Tr¥R(X + 5X) 
We will now derive an explicit relation between the generalized density 



(A3) 



R(X) = {0) R 
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and its perturbed value, expanded to first order in SX, 

R(X + SX) = V>R + V)R 
Note that the idempotence condition in Eq. ((5]) implies that the matrix Mr has the form 

m~ I (1) ^ 12 \ , ^ 

Mr = A4 

in the quasiparticle representation that diagonalizes (°>R. A straightforward linearization of 
Eq. ((A2| about <®R gives the relation 

Mr = SXM~ l F (A5) 

where M is the QRPA matrix, whose elements are given by second-order derivatives of the 



energy with respect to the generalized density matrix 
the vector notation 



28l |. and where we have introduced 




(A6) 



and similarly for Mr. Next, from Eqs. flAll) and (1A3|) . we deduce 

Sf ee f(X + 5X)-f(X) 



^ ■ Mr (A7) 



Combining this result with Eq. (1A5I) . we can express 5X in the form 



SX = - f ; ^ (A8) 



Equations (IA5I 1 and ( 1A8j ) are the basis for the iterative procedure described in the next 
section that is used to solve the HFB equation under constraint. 

In order to obtain a computationally efficient expression for the inverse QRPA matrix 
M~ l in Eq. (IA8|) . we adopt the so-called "cranking" approximation where the residual 
interaction between quasiparticles is neglected in the QRPA matrix. In this case, M -1 takes 
the block-diagonal form 



M 



-i 



[(e^ + Sv) 1 S liff S VT ] [0] 

[0] [(£^ + £u)~ 1 S tUT 5 VT ] 
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and therefore, 




(A9) 



with a corresponding expression for SX. 

2. Adjustment of the HFB solution in the case of one constraint 

In this section, we examine in greater detail steps [7] and [8] in the description of the HFB 
algorithm listed in section III Dl In this case, the constrained HFB equation is written 



where / is the expectation value of the constraint operator. The solution of the HFB equation 
then consists not only in determining R, but also the Lagrange multiplier A that satisfies the 
constraint. To solve this problem, we are led to an iterative procedure wherein the Lagrange 
multiplier is adjusted at each iteration. Consider the n th iteration, such that the generalized 
density matrix obtained in the previous iteration is R^ n ~^ with a corresponding Lagrange 
multiplier A (n_1) . The diagonalization of H (R^' 1 ^) - A^~% leads to a new generalized 
density which we will denote R^ n \ At this stage, the constraint is no longer necessarily 
satisfied and we calculate the deviation from the desired value 



[H (R) - AF, R] = 



with 




S = f — 



We correct the Lagrange multiplier using Eq. (1A8I) 



A (n) = A (n-1) + 
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and the generalized density using Eq. (|A5l) 



R (n) = R (n) + 5X M- l F 
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with 



S\ = A (n) - A (ri_1) 

We define the n th iteration with the self-consistent pair of and X^ n \ Note that the 
constraint is satisfied at each iteration. This iterative process generally converges, i.e. 

R (n) _^ £(n) _^ ^ 

A (n) -> A 

/(») _» / 

If the difference in constraint values is very large between successive iterations (as may be 
the case in the first few iterations), the convergence rate can be improved by calculating the 
generalized density matrix at the n th iteration according to 

flW = (I -a) (# (n) + 5AM- 1 /) + aR {n ^ ] 

with the associated Lagrange multiplier 

AW = (1 - a) (A^ 1 ) + 5X) + aX {n - 1] 

where the weight a tends to zero as the solution converges. With this prescription, the 
convergence of the generalized density and Lagrange multiplier are slowed down by the 
same amount. In other words, the desired value / for the constraint is approached in a 
gradual manner, so that at the n th iteration 



A (n) 



X (n) ) = f n) = (l-a)f + af {n - 



i) 



3. Adjustment of the HFB solution in the case of multiple constraints 

The results in the previous section can be readily generalized to an arbitrary number N 
of constraints. In this case, the HFB procedure minimizes the energy 



{A} 



N 



i=i 



{A} 



subject to the set of constraints 



{A} Fi {A}> = f u i = l,...,N 
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The generalized density matrix is now a function of N Lagrange multipliers, R({\}). We 
write 



R({X + 5X}) -R({\}) = W R 

'dX 



i=i 

A' 



(A10) 



i=l 



Clearly, is a variation where all the Lagrange multipliers are held fixed except for the 
one associated with Fj. Therefore, is given by Eq. (1A5I) with the substitutions 5X — ► <5Aj 
and F — > Fj. In the case of multiple constraints, Eq. (IA5I) is therefore replaced by 

AT 

Mr = ^SXiM^Fi (All) 

i=l 

Furthermore, using the generalization of Eq. (1A7I) to multiple constraints, 



Sfi 



{X + SX} Ft {A + 5A})-({A} Fi {A} 

If} ■ v>r 



and taking into account Eq. (1A11I) . we finally obtain 



5X = T- l 5f 



(A12) 



where the N x N matrix T is defined by 

Tim = 



'F ) 

1 m J 



(A13) 



Note that this matrix introduces correlations between all the constraints. We assume in 
our discussion that the inverse matrix T~ l exists, i.e., that the constraints are independent. 
Eqs. (1A11I ) and ( 1A12I ) then replace Eqs. (1A5I) and (IA8I) in the adjustment method described 
above. 



Appendix B: TRANSLATION IN A FINITE HARMONIC OSCILLATOR BASIS 



In this section, we give the explicit form for the expansion of a translated harmonic- 
oscillator function in a harmonic-oscillator basis. We begin with the generating function for 
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the Cartesian harmonic-oscillator function (Eq. (A.l) in 35]). 

^+2^/6-^/(26*) = J b ^Y^_ t k$ k{x . b) (B1) 

Letting x — > x + Ax on both sides of Eq. (IBll) after some simplification, the left-hand side 
(LHS) can be written as 



LHS = \/b^e- A ^ x+Ax ^ b 

2 m+n/2 ( Ax / & ) 



00 00 

V" 1 " -I \ r 1 h) 

m+n 



X 

m=0 n=0 



where we have used Eq. (IBll) to express the LHS in terms of harmonic-oscillator functions. 
Equating like powers of the arbitrary variable t between the LHS and right-hand side (RHS), 
we obtain 

<5> k (x + Ax;b) = e - A ^+ A */ 2 )/ fe2 

A 2 m / 2 Vk\(Ax/b) m 
m\y/{k-m)\ 
x$ k _ m {x;b) (B2) 

This is still a finite sum over harmonic-oscillator functions, however an overall exponential 
factor depending on x remains, and must be eliminated in order to obtain the expansion 
of $fe (x + Ax; b) on the harmonic-oscillator basis. Thus, in general, we need to derive an 
expansion for the expression 

e 2ax/b2 $i(x;b) (B3) 

where a = —Ax/2 and i = k — m in our case. Starting from the generating function in 
Eq. (IBll) . and multiplying both sides by the exponential factor in Eq. (1B3I) . the LHS of Eq. 
( jBlh becomes after some simplification 

LHS = ^e^f^^e^^t + ^^frb) 

1=0 ^ 

Expanding in powers of the arbitrary variable t, this takes the form 

00 00 / 



lhs = V^-^EEE 

1=0 p=0 q=0 



q ) p\VV- 
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Therefore, equating like powers of t between LHS and RHS, we obtain 



oo I 

J2 /f,2 



1=0 q=0 
/ry\ l+k-2q 
(j) *lfa1>) 



q / {k-q)\VT\ 



x 



Using this result in Eq. ( IB2I ). we obtain 

$ fc (x + Ax; b) = e - A - 2 /K) 

x£^(-^) (B4) 



where 



Q(£) = 2 (fc+/)/2 J-y£ {k+iy2 



x 



EE 



m! (A; — m — q)\ 

m=0 q=0 K ^' 




Note that the expansion of the translated harmonic-oscillator function requires in principle 
and infinite number of terms. In practice, these translations are performed in a finite-sized 
basis, and the truncation of the sum in Eq. (IB4I) to those shells within the basis can lead 
to the appearance of tails for translated nuclear densities expanded in a finite harmonic- 
oscillator basis, as shown in Figs. [Q and El 
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